1 Obtención y preparación de los datos

Cargamos los arhivos contenidos en “infoClin_all_dataset.RData” y “CRC_stageII.RData” con los contenidos clínicos y de CRC.

# Cargamos el archivo infoClin_all_dataset.RData 16/03/2024
dataClin <- load("data/infoClin_all_dataset.RData")

# Cargamos el archivo CRC_stageII.RData 19/03/2024 dataCRC <-
# load('data/CRC_stageII.RData')

# Cargamos el archivo exp_mat.RData 02/04/2024
dataCRC <- load("data/exp_mat.RData")

GSE Files

  • GSE39582 (Marisa et al. 2013) 421 vs 585
  • GSE14333 (Jorissen et al. 2009) 185 vs 290
  • GSE13294 (Jorissen et al. 2008) 121 vs 155
  • GSE17536 (Smith et al. 2010) 111 vs 177

Platform GPL570: [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array

Fusionamos los datos en un único dataframe al que le añadimos la columna grupo con los datos stage y msi.

Factorizamos las variables categóricas y ordenamos por ID tanto los datos clínicos como los de la matriz de expresión.

Creamos un DataFrame para asignar cada ID a su estudio clínico. A cada estudio le asignamos un color para poder identificarlo fácilmente en los estudios gráficos.

# 11/04/2023 Utilizamos el paquete 'dplyr'
library(dplyr)

# Función para preparar el data frame con el estudio por ID. También
# asignaremos un color a cada estudio
prepare_df <- function(df, study_name, color) {
    df %>%
        select(ID) %>%
        mutate(study = study_name, color = color)
}

# Lista de estudios con sus colores asociados
studies_info <- list(clx = "red", GSE39582 = "blue", GSE14333 = "green", GSE13294 = "orange",
    GSE17536 = "yellow", TCGA = "violet")

# Creamos una lista de de los estudios utilizando la función prepare_df
df_studies_list <- lapply(names(studies_info), function(x) {
    prepare_df(get(paste0(x, "_infoClin")), x, studies_info[[x]])
})

# Juntamos los df's de la lista se studies en uno solo y lo ordenamos por 'ID'
df_studies <- bind_rows(df_studies_list) %>%
    arrange(ID)

# Fusionamos y ordenamos todos los datos clínicos generando 'df_dataClin'
df_dataClin <- bind_rows(lapply(names(studies_info), function(x) get(paste0(x, "_infoClin")))) %>%
    arrange(ID)

# Cramos una nueva columna 'grupo' a los 'df_dataClin' con los datos de msi y
# stage
df_dataClin$grupo <- with(df_dataClin, paste(stage, msi_imputed, sep = "_"))

# Convertimos a factores
df_dataClin <- df_dataClin %>%
    mutate(across(c(msi_imputed, cms, stage, grupo), factor))

# Añadimos las columnas 'study' y 'color' a df_dataClin
df_dataClin <- left_join(df_dataClin, select(df_studies, ID, study, color), by = "ID")

Pasamos la matriz de expresión a data_frame y la ordenamos por “ID”

# Pasamos la matriz de expresión a data frame 26/03/2024
df_exmat <- as.data.frame(ex)

# Ordenamos la matriz de expresión por nombre de columnas 26/03/2024
df_exmat <- df_exmat[, order(names(df_exmat))]

2 Exploración de los datos

Resumen exploratorio de los datos:

# Hacemos un Summary con el paquete skimr 18/03/2024 install.packages('skimr')
require(skimr)
skim(df_dataClin)
Data summary
Name df_dataClin
Number of rows 1079
Number of columns 7
_______________________
Column type frequency:
character 3
factor 4
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
ID 0 1 7 12 0 1079 0
study 0 1 3 8 0 6 0
color 0 1 3 6 0 6 0

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
msi_imputed 0 1.00 FALSE 2 MSS: 840, MSI: 239
cms 158 0.85 FALSE 4 CMS: 329, CMS: 262, CMS: 194, CMS: 136
stage 0 1.00 FALSE 2 II: 684, III: 395
grupo 0 1.00 FALSE 4 MSS: 511, MSS: 329, MSI: 173, MSI: 66

Estadísticas y cálculos para la generación de los gráficos de datos clínicos:

# Estadísticas desglosadas 18/03/2024 Frecuencias
msi_frecuencias <- table(df_dataClin$msi_imputed)
cms_frecuencias <- table(df_dataClin$cms)
stage_frecuencias <- table(df_dataClin$stage)

# Valores NA
na_msi <- sum(is.na(df_dataClin$msi_imputed))
na_cms <- sum(is.na(df_dataClin$cms))
na_stage <- sum(is.na(df_dataClin$stage))

# Contamos cuantos MSI y MSS tenemos por cada estadio, II y III 19/03/2024
sII_MSI <- sum(df_dataClin$msi_imputed == "MSI" & df_dataClin$stage == "II")
sII_MSS <- sum(df_dataClin$msi_imputed == "MSS" & df_dataClin$stage == "II")
sIII_MSI <- sum(df_dataClin$msi_imputed == "MSI" & df_dataClin$stage == "III")
sIII_MSS <- sum(df_dataClin$msi_imputed == "MSS" & df_dataClin$stage == "III")

# Crear un dataframe con estos valores
stage_MSI_MSS_df <- data.frame(Stage_MSI_MSS = c("II_MSI", "II_MSS", "III_MSI", "III_MSS"),
    Freq = c(sII_MSI, sII_MSS, sIII_MSI, sIII_MSS))

3 Análisis de los datos clínicos

3.1 MSI/MSS

# Cargamos la librería "ggplot"
require(ggplot2)

# Pasamos a DataFrame msi
df_msi_frecuencias <- as.data.frame(msi_frecuencias)
df_msi_frecuencias
  Var1 Freq
1  MSI  239
2  MSS  840
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(df_msi_frecuencias$Freq)
# Añadimos una columna de porcentajes al dataframe
df_msi_frecuencias$Porcentaje <- (df_msi_frecuencias$Freq / total_frecuencias) * 100

# Diagrama de Barras de MSI/MSS
p1 <- ggplot(df_msi_frecuencias, aes(x = Var1, y = Freq)) +
  geom_bar(stat = "identity", fill = c("salmon", "seashell3")) +
  geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
  # Porcentajes en la mitad + o -
  geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 5, size = 3.5) +
  labs(title = "MSI/MSS", x = "", y = "Frecuencia") +
  theme_minimal()
p1

3.2 CMS

# Pasamos a DataFrame cms
df_cms_frecuencias <- as.data.frame(cms_frecuencias)
df_cms_frecuencias
  Var1 Freq
1 CMS1  194
2 CMS2  329
3 CMS3  136
4 CMS4  262
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(df_cms_frecuencias$Freq)
# Añadimos una columna de porcentajes al dataframe
df_cms_frecuencias$Porcentaje <- (df_cms_frecuencias$Freq / total_frecuencias) * 100

# Diagrama de Barras de CMS
p2 <- ggplot(df_cms_frecuencias, aes(x = Var1, y = Freq, fill = Var1)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = Freq), vjust = -0.3, size = 3.5) +  # Frecuencias arriba de cada barra
  # Porcentajes en la mitad + o -
  geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 5, size = 3.5) +  
  scale_fill_manual(values = c("salmon", "seashell3", "turquoise3", "slategray")) +
  labs(title = "CMS", x = "", y = "Frecuencia") +
  theme_minimal() +
  theme(legend.position = "none") # Omitir la leyenda

p2

3.3 ESTADIO

# Pasamos a DataFrame stage
df_stage_frecuencias <- as.data.frame(stage_frecuencias)
df_stage_frecuencias
  Var1 Freq
1   II  684
2  III  395
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(df_stage_frecuencias$Freq)
# Añadimos una columna de porcentajes al dataframe
df_stage_frecuencias$Porcentaje <- (df_stage_frecuencias$Freq / total_frecuencias) * 100

# Diagrama de Barras de ESTADIO
p3 <- ggplot(df_stage_frecuencias, aes(x = Var1, y = Freq)) +
  geom_bar(stat = "identity", fill = c("salmon", "seashell3")) +
  geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
  # Porcentajes en la mitad + o -
  geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 5, size = 3.5) +
  labs(title = "STAGE", x = "", y = "Frecuencia") +
  theme_minimal()
p3

3.4 ESTADIO vs MSI/MSS

# Stage vs MSI/MSS DataFrame
stage_MSI_MSS_df
  Stage_MSI_MSS Freq
1        II_MSI  173
2        II_MSS  511
3       III_MSI   66
4       III_MSS  329
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(stage_MSI_MSS_df$Freq)
# Añadimos una columna de porcentajes al dataframe
stage_MSI_MSS_df$Porcentaje <- (stage_MSI_MSS_df$Freq / total_frecuencias) * 100

# Diagrama de Barras de Stage vs MSI/MSS
p4 <- ggplot(stage_MSI_MSS_df, aes(x = Stage_MSI_MSS, y = Freq)) +
  geom_bar(stat = "identity", fill = c("salmon", "seashell3",
                                       "turquoise3", "slategray")) +
  geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
  # Porcentajes en la mitad + o -
  geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 2, size = 3.5) +
  labs(title = "STAGE vs MSI/MSS", x = "", y = "Frecuencia") +
  theme_minimal()
 
p4

3.5 TODOS

# Múltiple Layout con ggplot
require(patchwork)
# En un grid 2x2
p1 + p2 + p3 + p4 + plot_layout(ncol = 2, nrow = 2)

4 Análisis PCA

Para el Análisis de Componentes Principales (PCA), normalizamos y escalamos (estandarizamos) la matriz de expresión. En nuestro caso, los datos ya están normalizamos (en principio con \(log_2\)) por lo que tan solo los escalamos.
La estandarización se realiza restando la media de cada variable y dividiendo el resultado por la desviación estándar de esa variable.

\[ Z = \frac{x-\mu}{\sigma} \] Utilizamos la función scale para escalar. A continuación realizamos el PCA con la transpuesta de la matriz de expresión dado que queremos un análisis de las muestras en función de los genes y no al revés.
IMP: al hacer el escalado antes de la transposición estamos estandarizando por genes y no por muestras.

# Escalamos/normalizamos la matriz de expresión. También centra. 28/04/2024 Al
# hacer el escalado antes de la transposición estamos escalando por genes y no
# por muestras.
df_exmat_scaled <- scale(df_exmat)
# Hacemos el PCA. Hay que hacer la transpuesta de la matriz de expresión (las
# columnas representan genes y las filas representen muestras), así los
# componentes principales reflejan la variabilidad entre muestras. 28/04/2024
# Na haría falta center = TRUE ya que al escalar ya centralizamos.
pca_result <- prcomp(t(df_exmat_scaled), center = TRUE, scale. = FALSE)

RESUMEN PCA

# Obtenemos los valores propios (la varianza explicada por cada componente)
varianza <- pca_result$sdev^2

# Calculamos la proporción de varianza explicada por cada componente
prop_varianza <- varianza/sum(varianza)

# Calculamos la proporción acumulativa de varianza explicada
prop_acumulada <- cumsum(prop_varianza)

# Crear un dataframe con esta información para las primeras 10 PC
resumen_pca <- data.frame(SD = pca_result$sdev[1:10], Varianza = prop_varianza[1:10],
    VarAcumulada = prop_acumulada[1:10])

print(resumen_pca)
          SD   Varianza VarAcumulada
1  13.721545 0.08197744   0.08197744
2  11.124090 0.05387873   0.13585617
3  10.811433 0.05089263   0.18674880
4   9.600035 0.04012675   0.22687555
5   6.928531 0.02090117   0.24777672
6   6.605919 0.01900005   0.26677677
7   6.062676 0.01600358   0.28278035
8   5.642600 0.01386267   0.29664302
9   5.370890 0.01255975   0.30920276
10  5.178365 0.01167545   0.32087821

La baja varianza explicada por los primeros dos componentes, PC1 (8,2%) y PC2 (5.4%), sugiere que los datos tienen una distribución de varianza muy dispersa y que no hay unas pocas direcciones principales de variabilidad que dominen.
Esto es coherente con el hecho de que estamos tratando con datos biológicos, en concreto de CRC, que suelen ser de alta dimensionalidad con múltiples factores contribuyendo a la variabilidad general. Es decir, no hay patrones dominantes de variabilidad que puedan ser fácilmente capturados por los primeros componentes principales.

CARGAS DE GENES EN PCA

Para entender mejor qué genes están contribuyendo a la separación observada, podemos analizar las cargas de los componentes principales para ver qué genes tienen mayores pesos en los componentes que diferencian los grupos.

# Extracción de cargas de componentes principales
cargas_pca <- pca_result$rotation

# Mostramos las cargas de los primeros dos componentes
primer_componente <- cargas_pca[, 1]
segundo_componente <- cargas_pca[, 2]

PC1

# Extraemos los top 10 genes y sus cargas de PC1
top_genes_pc1 <- head(sort(primer_componente, decreasing = TRUE), 10)

# Convertimos a data frame para visualización de PC1
top_genes_pc1_df <- data.frame(Gen = names(top_genes_pc1), Carga = top_genes_pc1)
# Ordenamos el dataframe por carga de manera descendente
top_genes_pc1_df <- top_genes_pc1_df[order(top_genes_pc1_df$Carga, decreasing = TRUE),
    ]
# Ajustamos los factores para que el orden de los genes sea según las cargas
top_genes_pc1_df$Gen <- factor(top_genes_pc1_df$Gen, levels = top_genes_pc1_df$Gen)

# Creamos el gráfico de barras
plot_pc1 <- ggplot(top_genes_pc1_df, aes(x = Gen, y = Carga, fill = Gen)) + geom_bar(stat = "identity",
    position = "dodge") + theme_minimal() + theme(axis.text.x = element_text(angle = 90,
    hjust = 1, vjust = 0.5)) + labs(title = "Top 10 Cargas de Genes en PC1", x = "Gen",
    y = "Carga") + scale_color_gradient(low = "peachpuff", high = "peachpuff4")

# Mostramos los top 10 genes de PC1
print(top_genes_pc1)
     SFRP2       GAS1    COL10A1      CCL18     CYP1B1       SPP1      GPNMB 
0.06248246 0.04743630 0.04494118 0.04386073 0.04165178 0.03902618 0.03896303 
     GREM1     CCDC80      POSTN 
0.03860923 0.03848662 0.03844061

PC2

# Extraemos los top 10 genes y sus cargas de PC2
top_genes_pc2 <- head(sort(segundo_componente, decreasing = TRUE), 10)

# Convertimos a data frame para visualización de PC2
top_genes_pc2_df <- data.frame(Gen = names(top_genes_pc2), Carga = top_genes_pc2)
# Ordenamos el dataframe por carga de manera descendente
top_genes_pc2_df <- top_genes_pc2_df[order(top_genes_pc2_df$Carga, decreasing = TRUE),
    ]
# Ajustamos los factores para que el orden de los genes sea según las cargas
top_genes_pc2_df$Gen <- factor(top_genes_pc2_df$Gen, levels = top_genes_pc2_df$Gen)

# Creamos el gráfico de barras
plot_pc2 <- ggplot(top_genes_pc2_df, aes(x = Gen, y = Carga, fill = Gen)) + geom_bar(stat = "identity",
    position = "dodge") + theme_minimal() + theme(axis.text.x = element_text(angle = 90,
    hjust = 1, vjust = 0.5)) + labs(title = "Top 10 Cargas de Genes en PC2", x = "Gen",
    y = "Carga") + scale_color_gradient(low = "peachpuff", high = "peachpuff4")

# Mostramos los top 10 genes de PC1
print(top_genes_pc2)
      REG4       ZIC2      REG1A      TRIM7       CTSE       TCN1       AGR3 
0.07056229 0.05049736 0.04890627 0.04815804 0.04715623 0.04230700 0.04170686 
    RAB27B    SDR16C5    PLA2G2A 
0.04125286 0.04028539 0.03911398 
# Mostramos gráficamente las cargas del top 10 tanto de PC1 como de PC2
library(patchwork)
# En un grid 2x1
plot_pc1 + plot_pc2 + plot_layout(ncol = 2, nrow = 1)

GRÁFICAS PCA

Utilizamos el paquete factoextra para realizar las gráficas de los 2 primeros componentes, en concreto cla función fviz_pca_ind, que además nos permite añadir las elipses de confianza.

4.1 Estudios de origen

Entre los diferentes grupos de estudio, las elipses de confianza están centradas y solapadas lo que sugiere que, en el espacio de los primeros dos componentes principales, no hay una separación clara entre los grupos.

No es descartable la necesidad de examinar componentes adicionales o utilizar otras técnicas de análisis como métodos de Clustering.

# Utilizamos "factoextra" para hacer la gráfica de la PCA
library(factoextra)
# Hacemos la gráfica de los dos primeros componentes
fviz_pca_ind(pca_result,
             col.ind = df_dataClin$study, # Colores de los puntos basados en los estudios
             label = "none", 
             addEllipses = TRUE, # Añadir elipses de confianza
             ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
             legend.title = "Study",
             palette = "jco", # Paleta de colores, 'jco' es solo una opción
             title = "PCA Estudios")

4.2 MSI/MSS

Cuando representamos la PCA por inestabilidad de microsatélites vemos que las dos elipses, las de MSI y las de MSS, están claramente separadas lo que sugiere que ambos grupos tienen perfiles de expresión genética distintos en las dimensiones capturadas por los componentes principales (CP).
Como ya sabemos, estas diferencias son biológicamente relevantes; los componentes que muestran esta separación probablemente están capturando las variaciones genéticas ya conocidas de antemano.

# Hacemos la gráfica de los dos primeros componentes en relación a los grupos 
fviz_pca_ind(pca_result,
             col.ind = df_dataClin$msi_imputed, # Colores de los puntos basados en MSI/MSS
             label = "none", 
             addEllipses = TRUE, # Añadir elipses de confianza
             ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
             legend.title = "MSI/MSS",
             palette = "jco", # Paleta de colores, 'jco' es solo una opción
             title = "PCA MSI/MSS")

4.3 ESTADIO II/III

Entre los dos grupos de estadio, II y III, las elipses de confianza están más o menos centradas y algo solapadas lo que sugiere que, no hay una separación clara entre los grupos. Algo un tanto previsible ya que el estadio del cáncer depende del tiempo de progreso de la enfermedad.

# Hacemos la gráfica de los dos primeros componentes en relación a los grupos 
fviz_pca_ind(pca_result,
             col.ind = df_dataClin$stage, # Colores de los puntos basados en Estadio
             label = "none", 
             addEllipses = TRUE, # Añadir elipses de confianza
             ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
             legend.title = "Estadio",
             palette = "jco", # Paleta de colores, 'jco' es solo una opción
             title = "PCA Estadio II/III")

4.4 MSI/MSS vs ESTADIO II/III

Se podría decir que esta gráfica és una combinación de las dos anteriores; fijándonos en las elipses de confianza, se observan dos grupos separados claramente de otros dos, debido principalmente a la contribución de la instabilidad de microsatélites y no tanto a la del estadio del CRC.

# Hacemos la gráfica de los dos primeros componentes en relación a los grupos 
fviz_pca_ind(pca_result,
             col.ind = df_dataClin$grupo, # Colores de los puntos basados en Grupos
             label = "none", 
             addEllipses = TRUE, # Añadir elipses de confianza
             ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
             legend.title = "Grupo",
             palette = "jco", # Paleta de colores, 'jco' es solo una opción
             title = "PCA Grupos")

4.5 CMS

Al igual que con la inestabilidad de microsatélites, cuando representamos la PCA por clasificación CMS vemos que las cuatro elipses están claramente separadas; los 4 grupos tienen perfiles de expresión genética distintos en las dimensiones capturadas por los CP.
Estas diferencias son biológicamente relevantes dado que los componentes que muestran esta separación están capturando variaciones genéticas ya conocidas.

# Hacemos la gráfica de los dos primeros componentes en relación a los grupos 
fviz_pca_ind(pca_result,
             col.ind = df_dataClin$cms, # Colores de los puntos basados en CMS
             label = "none", 
             addEllipses = TRUE, # Añadir elipses de confianza
             ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
             legend.title = "CMS",
             palette = "jco", # Paleta de colores, 'jco' es solo una opción
             title = "PCA por CMS")

5 Estimación de la infiltración celular

Analizamos los datos viendo la infiltración celular.

5.1 MCPcounter

Utilizamos el paquete MCPcounter (Becht et al. 2016) que estima la abundancia poblacional de poblaciones de células inmunes (8 tipos) y estromales (2 tipos) que se infiltran en tejidos mediante la expresión genética.

En concreto, utilizatemos MCPcounter para estimar la infiltración celular del estroma y el sistema inmunitario en el microambiente tumoral (TME).

# Cargamos la librería 'MCPcounter'
require(MCPcounter)

# Con todos los datos de la matriz de expresión 26/03/2024
exmat_MCP <- MCPcounter.estimate(expression = df_exmat, featuresType = "HUGO_symbols")
# Transponemos
exmat_MCP <- t_df(exmat_MCP)
# Añadimos los datos clínicos a los datos de la matriz generada con MCPcounter
# 26/03/2024
exmat_MCP_Clin <- cbind(df_dataClin, exmat_MCP)

A modo de recordatorio didáctico de Inferencia Estadística:

  • El nivel de significación \(\alpha\) de un contraste es el error máximo de tipo I que estamos dispuestos a asumir.
\(H_0\) es verdadera \(H_0\) es falsa
Rechazar \(H_0\) Error Tipo I (α) Decisión Correcta
No rechazar \(H_0\) Decisión Correcta Error Tipo II (β)

  • El p-valor es la probabilidad del resultado del estadístico de contraste observado o de uno más alejado cuando la hipótesis nula es cierta.

  • El p-valor asociado a una observación del estadístico de contraste es el menor nivel de significación que nos permite rechazar la hipótesis nula.

  • Si el p-valor es inferior al nivel de significación \(\alpha\), rechazaremos la hipótesis nula, en caso contrario la aceptaremos.

En nuestro caso utilizamos un nivel de significación \(\alpha = 0.05\)

5.1.1 Infiltración MSI/MSS vs stage

5.1.1.1 MSI/MSS vs stage de todos los tipos celulares

Para hacer una comparativa entre los 4 grupos, utilizaremos el test no paramétrico de Kruskal-Wallis. Se suele utilizar para determinar si hay diferencias significativas entre tres o más grupos independientes cuando los datos no cumplen los supuestos necesarios para realizar un ANOVA: normalidad y homocedasticidad (homogeneidad de varianzas).

# Calculamos los test no paraméticos Kruskal-Wallis (asumimos que no tenemos
# normalidad) 16/04/2024

# tests de Kruskal-Wallis para NK
kw_NK_result <- kruskal.test(exmat_MCP_Clin$`NK cells` ~ grupo, data = exmat_MCP_Clin)

# Para el resto de tipos celulares utilizamos una lista con los nombres
tipos_celulares <- c("NK cells", "T cells", "CD8 T cells", "Cytotoxic lymphocytes",
    "B lineage", "Monocytic lineage", "Myeloid dendritic cells", "Neutrophils", "Endothelial cells",
    "Fibroblasts")

# Aplicamos el test de Kruskal-Wallis a cada tipo celular con 'lapply'
resultados_kw <- lapply(tipos_celulares, function(cell_type) {
    kruskal.test(reformulate("grupo", response = cell_type), data = exmat_MCP_Clin)
})

# Añadimos los nombres de los tipos celulares
names(resultados_kw) <- tipos_celulares

# Extraemos los valores de p de todos los tests realizados
p_valores <- sapply(resultados_kw, function(x) x$p.value)

La hipótesis nula del test de Kruskal-Wallis \(H_0\) establece que todas las poblaciones (grupos) tienen distribuciones idénticas, es decir, las medianas de todos los grupos son iguales.

\[ H_0: \forall i, j \, (i \neq j) \Rightarrow \mu_i = \mu_j \]

La hipótesis alternativa \(H_1\) establece que al menos uno de los grupos tiene una distribución diferente en cuanto a la mediana, comparada con las otras poblaciones sin especificar cuál de los grupos es diferente.

\[ H_1: \exists i \neq j \mid \mu_i \neq \mu_j \]

# Generamos todos los boxplot 27/03/2024
grupo <- exmat_MCP_Clin$grupo
# NK
bp1 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`NK cells`, grupo,
    "NK cells", p_valores["NK cells"])
# T cells
bp2 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`T cells`, grupo,
    "T cells", p_valores["T cells"])
# CD8 T cells
bp3 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`CD8 T cells`, grupo,
    "CD8 T cells", p_valores["CD8 T cells"])
# Cytotoxic lymphocytes
bp4 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Cytotoxic lymphocytes`,
    grupo, "Cytotoxic lymphocytes", p_valores["Cytotoxic lymphocytes"])
# B lineage
bp5 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`B lineage`, grupo,
    "B lineage", p_valores["B lineage"])
# Monocytic lineage
bp6 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Monocytic lineage`,
    grupo, "Monocytic lineage", p_valores["Monocytic lineage"])
# Myeloid dendritic cells
bp7 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Myeloid dendritic cells`,
    grupo, "Myeloid dendritic cells", p_valores["Myeloid dendritic cells"])
# Neutrophils
bp8 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$Neutrophils, grupo,
    "Neutrophils", p_valores["Neutrophils"])
# Endothelial cells
bp9 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Endothelial cells`,
    grupo, "Endothelial cells", p_valores["Endothelial cells"])
# Fibroblasts
bp10 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$Fibroblasts, grupo,
    "Fibroblasts", p_valores["Fibroblasts"])

library(patchwork)
# En un grid 3x2
bp1 + bp2 + bp3 + plot_layout(ncol = 3, nrow = 1)

bp4 + bp5 + bp6 + plot_layout(ncol = 3, nrow = 1)

bp7 + bp8 + bp9 + plot_layout(ncol = 3, nrow = 1)

bp10 + plot_layout(ncol = 3, nrow = 1)

5.1.1.2 Infiltración NKs

NORMALIDAD Y HOMOCEDASTICIDAD

Primero compromabos la normalidad de los grupos para las NK cells.

Tabla tabla con los resultados de los tests de Shapiro para normalidad.
Si p-value < 0.05 rechazamos la hipótesis nula de igualdad de estadísticos.

Shapiro test MSI MSS II III
p-valor 3.339e-05 1.058e-10 1.371e-10 1.089e-04

Dado que no hay normalidad no hace falta comprobar homocedasticidad.

Para realizar una comparativa entre dos grupos en los que asumimos no normalidad , en nuestro caso MSI vs MSS y estadio II vs estadio III, utilizamos el test no paramétrico de Mann-Whitney U también conocido como el test de Wilcoxon en el que se asume que los datos no necesitan ser normalmente distribuidos aunque deben poder ordenarse.

# Para dos grupos utilizamos el test no paramético Mann-Whitney U (asumimos que
# no tenemos normalidad) 16/04/2024

# Realizamos el test de Mann-Whitney U para los grupos MSI y MSS
mwU_NK_msi_test <- wilcox.test(grupo_msi, grupo_mss, alternative = "two.sided")

# Realizamos el test de Mann-Whitney U para los estadios II y III
mwU_NK_stage_test <- wilcox.test(grupo_sII, grupo_sIII, alternative = "two.sided")

La hipótesis nula \(H_0\) para el test de Mann-Whitney U establece que no hay diferencia en las medianas (o, más generalmente, en las distribuciones) de los dos grupos: \[ H_0: \mu_{grupo1} = \mu_{grupo2} \] La hipótesis alternativa \(H_1\) establece que hay una diferencia entre las distribuciones de los dos grupos. En nuestro caso la planteamos Bilateral (Two-sided):

\[ H_1: \mu_{grupo1} \neq \mu_{grupo2} \]

# Generamos todos los boxplot referentes a NKs 27/03/2024 Diagrama de Barras de
# NKs MSI vs MSS
bp11 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`NK cells`,
    exmat_MCP_Clin$msi_imputed, "NK cells MSI vs MSS", mwU_NK_msi_test$p.value)
# Diagrama de Barras de NKs II vs III
bp12 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`NK cells`,
    exmat_MCP_Clin$stage, "NK cells II vs III", mwU_NK_stage_test$p.value)

require(patchwork)
# En un grid 2x2
p4 + bp1 + p1 + bp11 + p3 + bp12 + plot_layout(ncol = 2, nrow = 3)

5.1.2 Infiltración resto de tipos celulares

5.1.2.1 T cells

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp13 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`T cells`,
    exmat_MCP_Clin$msi_imputed, "T cells MSI vs MSS")
# Diagrama de Barras II vs III
bp14 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`T cells`,
    exmat_MCP_Clin$stage, "T cells II vs III")

require(patchwork)
# En un grid 1x3
bp2 + bp13 + bp14 + plot_layout(ncol = 3, nrow = 1)

5.1.2.2 CD8 T cells

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp15 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`CD8 T cells`,
    exmat_MCP_Clin$msi_imputed, "CD8 T cells MSI vs MSS")
# Diagrama de Barras II vs III
bp16 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`CD8 T cells`,
    exmat_MCP_Clin$stage, "CD8 T cells II vs III")

require(patchwork)
# En un grid 1x3
bp3 + bp15 + bp16 + plot_layout(ncol = 3, nrow = 1)

5.1.2.3 Cytotoxic lymphocytes

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp17 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Cytotoxic lymphocytes`,
    exmat_MCP_Clin$msi_imputed, "Cytotoxic lymphocytes MSI vs MSS")
# Diagrama de Barras II vs III
bp18 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Cytotoxic lymphocytes`,
    exmat_MCP_Clin$stage, "Cytotoxic lymphocytes II vs III")

require(patchwork)
# En un grid 1x3
bp4 + bp17 + bp18 + plot_layout(ncol = 3, nrow = 1)

5.1.2.4 B lineage

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp19 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`B lineage`,
    exmat_MCP_Clin$msi_imputed, "B lineage MSI vs MSS")
# Diagrama de Barras II vs III
bp20 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`B lineage`,
    exmat_MCP_Clin$stage, "B lineage II vs III")

require(patchwork)
# En un grid 1x3
bp5 + bp19 + bp20 + plot_layout(ncol = 3, nrow = 1)

5.1.2.5 Monocytic lineage

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp21 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Monocytic lineage`,
    exmat_MCP_Clin$msi_imputed, "Monocytic lineage MSI vs MSS")
# Diagrama de Barras II vs III
bp22 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Monocytic lineage`,
    exmat_MCP_Clin$stage, "Monocytic lineage II vs III")

require(patchwork)
# En un grid 1x3
bp6 + bp21 + bp22 + plot_layout(ncol = 3, nrow = 1)

5.1.2.6 Myeloid dendritic cells

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp23 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Myeloid dendritic cells`,
    exmat_MCP_Clin$msi_imputed, "Myeloid dendritic cells MSI vs MSS")
# Diagrama de Barras II vs III
bp24 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Myeloid dendritic cells`,
    exmat_MCP_Clin$stage, "Myeloid dendritic cells II vs III")

require(patchwork)
# En un grid 1x3
bp7 + bp23 + bp24 + plot_layout(ncol = 3, nrow = 1)

5.1.2.7 Neutrophils

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp25 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$Neutrophils,
    exmat_MCP_Clin$msi_imputed, "Neutrophils MSI vs MSS")
# Diagrama de Barras II vs III
bp26 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$Neutrophils,
    exmat_MCP_Clin$stage, "Neutrophils II vs III")

require(patchwork)
# En un grid 1x3
bp8 + bp25 + bp26 + plot_layout(ncol = 3, nrow = 1)

5.1.2.8 Endothelial cells

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp27 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Endothelial cells`,
    exmat_MCP_Clin$msi_imputed, "Endothelial cells MSI vs MSS")
# Diagrama de Barras II vs III
bp28 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Endothelial cells`,
    exmat_MCP_Clin$stage, "Endothelial cells II vs III")

require(patchwork)
# En un grid 1x3
bp9 + bp27 + bp28 + plot_layout(ncol = 3, nrow = 1)

5.1.2.9 Fibroblasts

# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp29 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$Fibroblasts,
    exmat_MCP_Clin$msi_imputed, "Fibroblasts MSI vs MSS")
# Diagrama de Barras II vs III
bp30 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$Fibroblasts,
    exmat_MCP_Clin$stage, "Fibroblasts II vs III")

require(patchwork)
# En un grid 1x3
bp10 + bp29 + bp30 + plot_layout(ncol = 3, nrow = 1)

5.1.2.10 NK cells

# En un grid 1x3
bp1 + bp11 + bp12 + plot_layout(ncol = 3, nrow = 1)

5.1.3 Infiltración CMS

La clasificación CMS divide los tumores de CRC en cuatro subtipos principales, basados en características moleculares distintas, lo que refleja diferencias en la biología del tumor, la interacción con el TME, el pronóstico y la respuesta a tratamientos específicos. (Roelands et al. 2017)

Subtipo Características Inmunidad Prevalencia
CMS1 (Subtipo Inmuno) Alta inmunogenicidad, frecuentemente asociada con MSI y una alta carga mutacional. Tienen una infiltración activa de células inmunitarias efectoras, lo que puede hacerlos más susceptibles a terapias inmunogénicas. Aproximadamente el 14% de los casos de CRC.
CMS2 (Subtipo Canónico) Tumores con alta expresión de genes relacionados con el ciclo celular y la señalización WNT/MYC, presentando una marcada proliferación celular. Son generalmente poco inmunogénicos, con menor infiltración inmune comparado con CMS1 y CMS4. Es el subtipo más común, abarcando alrededor del 37% de los casos de CRC.
CMS3 (Subtipo Metabólico) Presenta alteraciones metabólicas, incluyendo un metabolismo energético distintivo, con una menor carga mutacional en comparación con CMS1. Similar a CMS2, muestra baja inmunogenicidad. Constituye aproximadamente el 13% de los casos de CRC
CMS4 (Subtipo Mesenquimal) Se caracteriza por la activación de vías relacionadas con la angiogénesis, la activación de células estromales y la señalización TGF-β, lo que conduce a un entorno inmunosupresor. Aunque presenta una significativa infiltración inmune, esta es predominantemente de naturaleza supresora, con una alta presencia de células estromales. Aproximadamente el 23% de los casos de CRC caen en este subtipo.

Vemos la infiltración celular según CMS.

# Calculamos los test no paraméticos Kruskal-Wallis (asumimos que no tenemos
# normalidad) 16/04/2024

# Aplicamos el test de Kruskal-Wallis a cada tipo celular con 'lapply'
resultados_kw_cms <- lapply(tipos_celulares, function(cell_type) {
    kruskal.test(reformulate("cms", response = cell_type), data = exmat_MCP_Clin)
})

# Añadimos los nombres de los tipos celulares
names(resultados_kw_cms) <- tipos_celulares

# Extraemos los valores de p de todos los tests realizados
p_valores_cms <- sapply(resultados_kw_cms, function(x) x$p.value)
# Quitamos los NA's de la columna cms
exmat_MCP_Clin_CMS <- exmat_MCP_Clin[!is.na(exmat_MCP_Clin$cms), ]
# Generamos todos los boxplot 02/04/2024
cms <- exmat_MCP_Clin_CMS$cms
# NK
bp1 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`NK cells`,
    cms, "NK cells", p_valores_cms["NK cells"])
# T cells
bp2 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`T cells`,
    cms, "T cells", p_valores_cms["T cells"])
# CD8 T cells
bp3 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`CD8 T cells`,
    cms, "CD8 T cells", p_valores_cms["CD8 T cells"])
# Cytotoxic lymphocytes
bp4 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Cytotoxic lymphocytes`,
    cms, "Cytotoxic lymphocytes", p_valores_cms["Cytotoxic lymphocytes"])
# B lineage
bp5 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`B lineage`,
    cms, "B lineage", p_valores_cms["B lineage"])
# Monocytic lineage
bp6 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Monocytic lineage`,
    cms, "Monocytic lineage", p_valores_cms["Monocytic lineage"])
# Myeloid dendritic cells
bp7 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Myeloid dendritic cells`,
    cms, "Myeloid dendritic cells", p_valores_cms["Myeloid dendritic cells"])
# Neutrophils
bp8 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$Neutrophils,
    cms, "Neutrophils", p_valores_cms["Neutrophils"])
# Endothelial cells
bp9 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Endothelial cells`,
    cms, "Endothelial cells", p_valores_cms["Endothelial cells"])
# Fibroblasts
bp10 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$Fibroblasts,
    cms, "Fibroblasts", p_valores_cms["Fibroblasts"])

require(patchwork)
# En un grid 3x2
p2 + bp1 + bp2 + plot_layout(ncol = 3, nrow = 1)

bp3 + bp4 + bp5 + plot_layout(ncol = 3, nrow = 1)

bp6 + bp7 + bp8 + plot_layout(ncol = 3, nrow = 1)

bp9 + bp10 + plot_layout(ncol = 3, nrow = 1)

Los tumores CMS2 y CMS3 son poco inmunogénicos; CMS1 y CMS4 difieren en el tipo de infiltración inmune.
Los tumores CMS1 tienden a acumular una gran cantidad de mutaciones debido al estado de MSI y atraen células efectoras inmunes.
Los tumores CMS4 exhiben un TME inmunosupresor con infiltración estromal. (Lanuza et al. 2022)

6 Referencias

Repositorio: (Ferrón 2024)

Becht, Etienne, Nicolas A Giraldo, Laetitia Lacroix, Bénédicte Buttard, Nabila Elarouci, Florent Petitprez, Janick Selves, et al. 2016. “Estimating the Population Abundance of Tissue-Infiltrating Immune and Stromal Cell Populations Using Gene Expression.” Genome Biology 17: 1–20.
Ferrón, Sheddad Kaid-Salah. 2024. “Caracterización Inmunonómica Del Microambiente Tumoral En Cáncer Colorrectal En Estadios Tempranos y Tardı́os.” https://github.com/Sheddad/TFM_Inmunonomica_TME_CRC.
Jorissen, Robert N, Peter Gibbs, Michael Christie, Saurabh Prakash, Lara Lipton, Jayesh Desai, David Kerr, et al. 2009. “Metastasis-Associated Gene Expression Changes Predict Poor Outcomes in Patients with Dukes Stage b and c Colorectal Cancer.” Clinical Cancer Research 15 (24): 7642–51.
Jorissen, Robert N, Lara Lipton, Peter Gibbs, Matthew Chapman, Jayesh Desai, Ian T Jones, Timothy J Yeatman, et al. 2008. “DNA Copy-Number Alterations Underlie Gene Expression Differences Between Microsatellite Stable and Unstable Colorectal Cancers.” Clinical Cancer Research 14 (24): 8061–69.
Lanuza, Pilar M, M Henar Alonso, Iratxe Uranga-Murillo, Sandra Garcı́a-Mulero, Cecilia Pesini, Eva M Gálvez, Ariel Ramı́rez-Labrada, Rebeca Sanz-Pamplona, and Julián Pardo. 2022. “Adoptive NK Cell Transfer as a Treatment in Colorectal Cancer Patients: Analyses of Tumour Cell Determinants Correlating with Efficacy in Vitro and in Vivo.” Frontiers in Immunology 13: 890836.
Marisa, Laetitia, Aurélien de Reyniès, Alex Duval, Janick Selves, Marie Pierre Gaub, Laure Vescovo, Marie-Christine Etienne-Grimaldi, et al. 2013. “Gene Expression Classification of Colon Cancer into Molecular Subtypes: Characterization, Validation, and Prognostic Value.” PLoS Medicine 10 (5): e1001453.
Roelands, Jessica, Peter JK Kuppen, Louis Vermeulen, Cristina Maccalli, Julie Decock, Ena Wang, Francesco M Marincola, Davide Bedognetti, and Wouter Hendrickx. 2017. “Immunogenomic Classification of Colorectal Cancer and Therapeutic Implications.” International Journal of Molecular Sciences 18 (10): 2229.
Smith, J Joshua, Natasha G Deane, Fei Wu, Nipun B Merchant, Bing Zhang, Aixiang Jiang, Pengcheng Lu, et al. 2010. “Experimentally Derived Metastasis Gene Expression Profile Predicts Recurrence and Death in Patients with Colon Cancer.” Gastroenterology 138 (3): 958–68.